if !d.name eq 'PS' then begin
    device,xsize=18,ysize=12,yoffset=3
    !p.thick=2 & !x.thick=1 & !y.thick=1
end
!p.multi=[0,2,1] & !p.charsize=1.3
!x.margin=[7,0] & !y.margin=[5,3]

@eu_setup
@rhoprof
;restore,'urms.sav'
;nn=where(r ge 0.72)
;urms=mean(surms[nn])
; density heigth scale
;hd=-1./xder_curl(alog(rh),r)
;tau=Hd/urms

;etat=10^4*rr*tau*surms^2/3.
ns=205D
ng=50D

rsinth1=fltarr(m,l)
dwdphi=fltarr(n,m,l,ns)
dvdphi=fltarr(n,m,l,ns)
dsinudth=fltarr(n,m,l,ns)
dwdth=fltarr(n,m,l,ns)
drvdr=fltarr(n,m,l,ns)
drudr=fltarr(n,m,l,ns)

wr=fltarr(n,m,l,ns)
wth=fltarr(n,m,l,ns)
wphi=fltarr(n,m,l,ns)

;TOTAL HELICITY COMPONENT

for it=0,ns-1 do begin
print,'Computing total helicity',it
 for k=0, l-1 do begin
  for i=0,n-1 do begin
   dsinudth[i,*,k,it]=xder_curl(reform(smooth(sin(theta)*u[i,*,k,it],1)),theta)       
   dwdth[i,*,k,it]=xder_curl(reform(smooth(w[i,*,k,it],1)),theta)
  endfor
 endfor

 for j=0,m-1 do begin
  for i=0,n-1 do begin
   drvdr[i,j,*,it]=xder_curl(smooth(r[*]*v[i,j,*,it],1),r)   
   drudr[i,j,*,it]=xder_curl(smooth(r[*]*u[i,j,*,it],1),r)
  endfor
 endfor
;

 for k=0, l-1 do begin
  for j=1,m-1 do begin
   rsinth1[j,k]=1./(r[k]*sin(theta[j]))
  endfor
 endfor
 rsinth1[0,*]=0.
 rsinth1[m-1,*]=0.


 for k=0, l-1 do begin
  for j=1,m-1 do begin
   for i=0,n-1 do begin
     wr[i,j,k,it]=rsinth1[j,k]*dsinudth[i,j,k,it]
     wth[i,j,k,it]=rsinth1[j,k]*dwdphi[i,j,k,it]-drudr[i,j,k,it]/r[k]
     wphi[i,j,k,it]=drvdr[i,j,k,it]/r[k]-dwdth[i,j,k,it]/r[k]
   endfor
  endfor
 endfor

endfor

Ht=total((wr*w + wth*v + wphi*u),1)/double(n)
Htm=total(Ht[*,*,ns-ng:ns-1],3)/double(ng)

; U x w, total
uxwt_r=v*wphi-u*wth
uxwt_th=u*wr-w*wphi
uxwt_phi=w*wth-v*wr
d_ruxwtth_dr=fltarr(n,m,l,ns)
d_uxwtr_dth=fltarr(n,m,l,ns)
curl_uxwt_phi=fltarr(n,m,l,ns)

; PHI COMPONENT OF curl U x w
for it=0,ns-1 do begin
 print,'computing phi component of Uxw total', it
 for j=0,m-1 do begin
  for i=0,n-1 do begin
   d_ruxwtth_dr[i,j,*,it]=xder_curl(smooth(r[*]*uxwt_th[i,j,*,it],1),r)   
  endfor
 endfor

 for k=0, l-1 do begin
  for i=0,n-1 do begin
   d_uxwtr_dth[i,*,k,it]=xder_curl(reform(smooth(uxwt_r[i,*,k,it],1)),theta)
  endfor
 endfor

 for k=0, l-1 do begin
  for j=1,m-1 do begin
   for i=0,n-1 do begin
    curl_uxwt_phi[i,j,k,it]=d_ruxwtth_dr[i,j,k,it]/r[k]-d_uxwtr_dth[i,j,k,it]/r[k]
   endfor
  endfor
 endfor

endfor
curl_uxwt_phi = total(curl_uxwt_phi,1)/double(n)
av_curl_uxwt_phi = total(curl_uxwt_phi[*,*,ns-ng:ns-1],3)/double(ng)
help,av_curl_uxwt_phi

; CURL OF THE MEAN
print,'Computing mean helicity',ns
uum=total(u,1)/double(n)
vvm=total(v,1)/double(n)
wwm=total(w,1)/double(n)

um=total(uum[*,*,ns-ng:ns-1],3)/double(ng)
vm=total(vvm[*,*,ns-ng:ns-1],3)/double(ng)
wm=total(wwm[*,*,ns-ng:ns-1],3)/double(ng)

dwdphim=fltarr(m,l)
dvdphim=fltarr(m,l)
dsinudthm=fltarr(m,l)
dwdthm=fltarr(m,l)
drvdrm=fltarr(m,l)
drudrm=fltarr(m,l)

wrm=fltarr(m,l)
wthm=fltarr(m,l)
wphim=fltarr(m,l)

 for k=0, l-1 do begin
   dsinudthm[*,k]=xder_curl(reform(smooth(sin(theta)*um[*,k],1)),theta)       
   dwdthm[*,k]=xder_curl(reform(smooth(wm[*,k],1)),theta)
 endfor
 for j=1,m-1 do begin
   drvdrm[j,*]=xder_curl(smooth(r[*]*vm[j,*],1),r)   
   drudrm[j,*]=xder_curl(smooth(r[*]*um[j,*],1),r)
 endfor
 for k=0, l-1 do begin
  for j=1,m-1 do begin
     wrm[j,k]=rsinth1[j,k]*dsinudthm[j,k]
     wthm[j,k]=rsinth1[j,k]*dwdphim[j,k]-drudrm[j,k]/r[k]
;     wthm[j,k]=-drudrm[j,k]/r[k]
     wphim[j,k]=drvdrm[j,k]/r[k]-dwdthm[j,k]/r[k]
  endfor
 endfor
; mean kinetic helicity
hm=wrm*wm + wthm*vm + wphim*um
; mean <Uxw>
uxwm_r=vm*wphim-um*wthm
uxwm_th=um*wrm-wm*wphim
uxwm_phi=wm*wthm-vm*wrm
d_ruxwmth_dr=fltarr(m,l)
d_uxwmr_dth=fltarr(m,l)
curl_uxwm_phi=fltarr(m,l)

; phi component of curl <Uxw>
 for j=1,m-1 do begin
   d_ruxwmth_dr[j,*]=xder_curl(smooth(r[*]*uxwm_th[j,*],1),r)   
 endfor
 for k=0, l-1 do begin
   d_uxwmr_dth[*,k]=xder_curl(reform(smooth(uxwm_r[*,k],1)),theta)
 endfor
 for k=0, l-1 do begin
  for j=1,m-1 do begin
   curl_uxwm_phi[j,k]=d_ruxwmth_dr[j,k]/r[k]-d_uxwmr_dth[j,k]/r[k]
  endfor
 endfor

; small scale quantities
khel=htm-hm

tkhel=transpose(khel)
rgood=0.95
dd=where(r le rgood)
thg=where(th GE -85 AND th LT 0)

curl_uxws = transpose(AV_CURL_UXWT_PHI-CURL_UXWM_PHI)
lev=grange(min(curl_uxws[thg,3:l-4]),max(curl_uxws[thg,3:l-4]),64)
contour,smooth(curl_uxws[*,thg],2),x[*,thg],y[*,thg],/fi,nl=64,/iso,lev=lev


save,filename='khelicity.sav',r,th,x,y,khel
;!p.multi=[0,2,1] & !p.charsize=1.3 
;contour,transpose(smooth(khel,3)),x,y,/fi,nl=64,/iso,title='!4a!8=-!4s!8<!4x!8`.u`>/3!6'
;plot,r,smooth(etat,3),/ylog,xtitle='!8r!6',ytitle='!4g!8!Dt!N!6(cm!U!82!N!6s!U!8-1!N!6)!6',$
;       title='!4g!8!Dt!N=!4s!8u!D!6rms!8!U2!N/3!6',thick=3


end
